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ABSTRACT 



With the advent of long term low Earth-oroiting satellites 
comes the requirement to maintain a specified orbital radial 
band. Optimal control theory suggests that periodic thrusting 
is more efficient than forced Keplerian motion in obital 
maintenance through the use of Primer vectoring. This thesis 
examined the efficiency of fixed-angle transverse thrusting as 
an alternative to Primer vectoring. Numerical analysis shows 
that a thrust angle of 65-70 degrees is feasible for radial 
band control for a wide range of parameters. Fuel usage can 
be minimized through the proper selection of radial bandwidth, 
thruster size, and thrust angle. This thesis shows that 
forced Keplerian motion is always superior to fixed-angle 
transverse thrusting from a fuel usage standpoint, and hence 
that thrust vectoring must be utilized. 
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INTRODUCTION 



Orbital maintenance for low-Earth orbiting (LEO) 
satellites can be accomplished using forced Keplerian motion 
or by bang-bang control. Forced Keplerian motion utilizes 
constant thrust, with magnitude equal to drag, to prevent 
orbital decay and maintain orbital altitude. Bang-bang 
control allows the satellite orbit to decay and then fires 
thrusters to reboost the satellite and keep it within a 
prescribed radial band. This can be modeled as a series of 
orbital transfers offsetting decay due to atmospheric drag. 
Currently, LEO satellites are not required to maintain a 
specified orbital band, and are reboosted only to prevent 
reentry or to change orbit. 

Historically, co-planar orbital transfers have been 
optimally accomplished using two or three-impulse Hohmann 
maneuvers [Ref. l:pp 78-88]. These maneuvers are designed 
for one-time transfers between exo-atmospher ic orbits, not 
for a series of maneuvers to offset drag and maintain an 
orbital band. Bang-bang control utilizes one thruster 
firing per reboost of the satellite, vice two or three 
firings used in a Hohmann transfer. 

In Ref. 2, the problem of minimum-fuel orbital 
maintenance for low-altitude, non-lifting bodies was 
considered, with the conclusion that fixed-angle, transverse 
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thrusting (i.e. bang-bang control with the thrust vector at 
a constant angle relative to the local horizontal) is less 
efficient than forced Keplerian motion. This topic was 
first explored by Ross and Melton [Ref. 3] who concluded 
that forced Keplerian motion is not the optimal solution for 
this problem. Development of long-term, low-altitude 
satellites, such as the space station, dictate continued 
research into this topic to resolve this apparent 
contradiction . 

Little research has been done on the optimization of 
bang-bang control in relation to atmospheric effects on non- 
lifting (blunt) bodies. A better understanding of these 
effects is essential for optimization of orbits and 
propulsion eguipment for long-term satellites, including the 
space station. Additionally, design and operation of 
lifting bodies may benefit from this study. 
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II. GENERAL FORMULATION OF THE PROBLEM 



Two methods of optimizing orbital transfers for LEO 
satellites have been studied in Ref . 2 and Ref. 3. On the 
basis of numerical simulations, Pauls [Ref. 2] concludes 
that bang-bang control is less efficient than forced 
Keplerian motion, and that transverse thrusting at a fixed 
angle of 70 degrees is generally the optimal feasible 
solution for orbital maintenance if bang-bang control is 
used. This apparently contradicts the analytical results of 
Ross and Melton [Ref. 3], who conclude that forced Keplerian 
motion is not optimal and that a bang-bang solution exists 
that optimizes fuel usage. This leads Pauls to conclude 
that variable thrust-vectoring is essential for optimality 
during finite-burns. However, Pauls study was limited to a 
small class of satellites (essentially space platforms) in 
the upper atmosphere. Therefore, his conclusions are not 
globally valid, and questions remain. Is transverse 
thrusting less efficient than forced Keplerian motion for 
all parameters? What exactly are these parameters? Is the 
70 degree thrust angle always optimal? This thesis will try 
to answer some of these questions. 

The problem is defined as maintaining a spacecraft 
within a prescribed radial band R,^ R^R^ . Comparison of 
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fuel usage for bang-bang control and forced Keplenan 
motion, over a wide range of parameters, will provide 
insights on the optimality question. Also, determination of 
a thrust angle, or range of angles, that optimizes the bang- 
bang control solution is desired. 

Ross and Melton [Ref. 3] base their conclusions on 
methods utilizing optimal control theory, which require 
initial and boundary conditions to facilitate a solution 
using numerical integration. This is a cumbersome solution 
technique that requires recomputation for slight changes in 
initial or boundary conditions. While this may be the most 
accurate solution technique, it is not easy or very 
practical . 

Pauls [Ref. 2] develops a computer model to solve this 
problem. Equations of motion were derived and solved for 
various parameter combinations to optimize bang-bang 
control. Attempts to define a radial band and then control 
the spacecraft within that band were unsuccessful. However, 
maintaining specific energy allowed a band to be maintained, 
albeit larger than the desired band. All solutions to this 
model predicted that forced Keplerian motion is more 
efficient than bang-bang control. 

Further study using a computer model is dictated. Non- 
dimensionalization of the equations of motion will allow 
effects of variations in parameters to be more easily 
analyzed. Variations of parameters will be analyzed for 
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trends and effects on solutions, leading to optimal 
solutions for design of bang-bang control systems ( sucn as 
thruster size, thrust angle, and size of radial band). 
Optimality will be based on fuel usage, maintenance of the 
reauirea radial band (with minimal eccentricity) , and 
thruster firing patterns. An ability to predict the actual 
radial band controlled is desired. 

Results of the computer model will then be compared to 
those obtained using forced Keplerian motion to determine 
the optimal means of maintaining a spacecraft within a 
design radial band while subjected to continuous atmospheric 
drag . 
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III. FORMULA DEVELOPMENT 



A. DEVELOPMENT OF THE EQUATIONS OF MOTION 

For simplification, orbital motion is assumed to be 
coplanar and initially circular. Assuming a non-lifting 
(blunt) body, drag becomes the only aerodynamic force 
perturbing the orbit and thrust is the external force 
applied to counter-balance drag. The external forces acting 
on the spacecraft consist of gravity, aerodynamic drag, and 
thrust. Figure 1 shows the coordinate system and net forces 
for this two-body problem. 



Transverse axis 




Earth 



Figure 1 Graphical Representation of Coordinate System 
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The equations of motion for this two-boav problem are 



m 



(l) 




m 



v/here a r and a^ are the radial and transverse accelerations 
and T F r and 2 F tr are the sums of the external radial and 
transverse forces, respectively, and m is spacecraft mass. 
From Figure 1, the components of drag are 



D r =-Dsin 
D, r =-D cos 

and the components of thrust are 



(7) 


(3) 


(7) 


(4) 



T r =Tsin(a) (5) 

T^=Tcos(a) (6) 

The angle y is the flight path angle and is defined by the 
intersection of the velocity vector and transverse axis, 
while angle a is the thrust angle and is defined by the 
intersection of the thrust vector and transverse axis. The 
equations of motion can now be written 



r 



n 



-9' 2 r= 




(7) 



d"r+2d'r' 




( 8 ) 
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where n is the Earth's gravitational constant, 9 is the 
angular position of the spacecraft, and primes denote 
differentiation with respect to time. 

B. NONDIMENS I ONALIZ ATI ON OF THE EQUATIONS OF MOTION 

The equations of morion are nondimensionalized so that 
variations in parameters are minimalized. The 
nondimensionalizing constants are chosen to "balance" the 
equations and thus produce the same order of magnitude 
responses. This facilitates study of the effects of 
changing parameters on the optimization problem. 

1. Definitions 

In order to nondimensionalize the equations, 
variables must be defined and nondimensionalizing constants 
must be chosen. Nondimensionalizing constants for radius 
and time are designated R and r, respectively, and were 
chosen somewhat arbitrarily based on the Earth's radius. 
Thus, R is the radius of the Earth and r is the period of a 
circular orbit whose radius equals R, and is defined as 

r =2n 

These constants define the nondimensionalized variables for 
radius and time 



& ( 9 ) 
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The nonaimensionalizing factor for mass is M (spacecraft 
initial mass), which defines the nondimensionalized mass 




The variable 6 is already dimensionless. 

2. Nondimensionalization 

Using the above definitions, the following 
relationships are developed. 

dr = t dr 

dZ ~R ~dt 

d 2 r _ d , r dr . _ r 2 d 2 r 
~ddd} ~di K ~R~dt‘ ~R~dt 2 

dd_ _ T de_ 

dZ dt 

£l = J_ (T d^ )=T 2d 2 (? 
dZ 2 dZ dt dt 2 

Substituting these equations into the equations 7 and 
substituting equation 9 for r where convenient, and 



( 12 ) 

(13) 

(14) 

(15) 

(16) 



rearranging terms, yields the following nondimensionalized 
results 




(17) 



d 2 d __ 2 dddr _ r 2 Dcosy + r 2 T cosa 
df 2 r df 2 Rr mM Rr mM 



(18) 



C. NONDIMENS IONALIZATION OF ATMOSPHERIC DRAG 



1. Definitions 



Drag is given by 




(19) 



where p is atmospheric density (altitude dependent) , S rcf is 
the reference surface area of the spacecraft affected by 
atmospheric density, C d is the coefficient of drag, and v is 
spacecraft velocity. 

2. Nondimensionalization 



Nondimensionalization of v 2 is accomplished using 



the relationship 



v 2 =r' 2 +r 2 6 ,2 = ( ) 2 +r 2 ( ^ ) 



( 20 ) 



Substituting equations 10,13, and 15 gives 
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which is then substituted into equation 19 and simplified 
into the form used in equations 17 and 13 to give 



T 2 D _ T 2 P S rer C d R ' V " _ pi?V 2 

R MW R 2MWt 2 2MB 



( 22 ) 



where B is the initial ballistic coefficient and defined as 



B=- 



M 



n c 



(23) 



This result is then put back into equations 17 and 18 to 
yield the nondimensionalized equations 



d 2 r J- 47r - pRv sin 7 + r 2 Tsina 

di 2 dX "r 2 " 2MB "R" MW 



(24) 



d 2 8 __ 2 d6dr _ pRv 2 cosy + r 2 Tcosa 
dT 2 r cffcff 2MBr i?r MW 



(25) 



Inspection of these equations reveals that further 
simplification is possible. Assuming an exponential 
atmospheric density model (which will be used throughout the 
study of this problem) , with atmospheric scale height jS 
defines 



P = Po e 






( 26 ) 



This allows the ballistic coefficient, B, to be 
nondimensionalized as 



( 27 ) 



Po* 

The thrust terms can be similarly nondimensionalized by 
defining a nondimensionalized thrust 




( 28 ) 



Substituting equations 27 and 28 into equations 24 and 25, 
and nondimensionalizing the exponent for drag result in the 
final equations 



d 2 r = - , dQ ' 2 . 4 tv 
dt 2 ~ dt r 2 



_e H v“smy 

2 mB 



Tsina 

m 



( 29 ) 



d 2 6 __ 2 dQdr 
dt 2 r dtdt 



-p«( ?--%£)_, 

_e ^ v~cosy + Tcosa 



( 30 ) 



2miB 



mi 



These equations will be used to generate the computer model 
to study the optimality problem of bang-bang control. The 
only parameters that need to be varied to study the 
optimality problem are 



a , P , T, B 

Since the only variable parameter in the nondimensionalized 
thrust term is thrust, T, either one can be varied for the 
study. 
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D. DEVELOPMENT OF MASS EQUATIONS 

One of the optimality criteria to be used in this study 
is che mass of fuel consumed to maintain orbit. The 
relationship 




(31) 



defines the rate of fuel consumption and can be integrated 
over time to determine the total fuel consumed. Using 
equations 11, 12, and 29, this can be nondimensionalized and 
simplified to 



dt x I sp g 

These equations will be used to calculate the fuel used to 
maintain orbit using bang-bang control for comparison to 
fuel used in forced Keplerian motion. 



13 



IV. DEVELOPMENT AND VALIDATION OF COMPUTER MODEL 



A. COMPUTER PROGRAM DEVELOPMENT 

A computer program, written in FORTRAN and listed in 
Appendix A with sample input and output files, was developed 
to simulate spacecraft orbital motion. A fourth-order 
Runge-Kutta numerical integration routine is used to 
integrate the equations of motion. The program is comprised 
of a main program section, which controls input, output, and 
flow of information, and five subroutines that provide 
computations necessary to simulate orbital motion and 
calculate the osculating orbital parameters. 

The non-dimensionalized equations previously developed 
(equations 29 and 30) are simplified for use in the computer 
program using the following state variable definitions. 



x,=r 


(33) 


d x \ _ dr 


(34) 


dt dt 




x 3 =0 


(35) 


dx 2 = d 0 


(36) 


dt dt 





These state variables are nondimensionalized as follows 



_ r *1 - 

X - = i = r- 

' R R 



( 37 ) 



_ = ar = r dr = r 
‘ 2 dZ R at R‘" 2 



(38) 



*4 



dZ dt 



= TX. 



(39) 



The state variable x 3 is already nondimensionalized. These 
equations are then substituted into equations 29 and 30 to 
give the equations of motion used in the computer program. 



dx2 

~dZ 



=x,x 4 ‘ - . 



47T 



" vsiny Tsina 



(40) 



X, 



2mE 



m 



2 — — e " rcosy Tcosq (41) 

2 — — — — 

Xj 2mrH mr 

The first subroutine, DRAG, calculates the atmospheric 
drag experienced by the spacecraft. A constant atmospheric 
density model is used for initial program validation, after 
which an exponential density model is used. The second 
subroutine, EQN, updates the equations of motion that define 
the spacecraft's orbital motion, equations 40 and 41. The 
third subroutine, ORBPAR, calculates the osculating orbital 
parameters for the spacecraft's motion, including semi-major 
axis, apogee, perigee, and period. The fourth subroutine, 
RK4, is a standard fourth-order Runge-Kutta numerical 



dxl 

~dZ 
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integration routine that integrates the state variables. 

The last subroutine, THRUST, governs the activation and 
deactivation of the thrusters to maintain the satellite's 
orbit . 

B. PROGRAM VALIDATION 

1. Initial validation with no external forces 

Program validation proceeded in a logical series of 
steps, with each step requiring positive validation before 
proceeding to the next step. Initially, all external forces 
except gravity were neglected in order to maintain a 
circular orbit. As shown in Appendix B, Figures 1-4, 
radius, velocity, angular momentum, and specific energy all 
remained constant. Initial conditions of an elliptic orbit 
were then input into the program with all external forces 
except gravity neglected. Figures 5-8 of Appendix B show 
that semi-major axis (SMA) , eccentricity, specific energy 
and angular momentum are constant. This validates the 
initial program. 

2. Validation of drag 

A constant atmospheric density was then introduced 
into the program, with a value that would produce a 
noticeable effect in a small number of orbits to enhance the 
validation process. The program was run for ten orbits and 
the results compared with analytical results obtained using 
equations from Ref. 4, 
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( 42 ) 



Av=-«r(_iJ*') pa v 



( 43 ) 



AP = -6n : (—^)p — 
M v 



( 44 ) 



which show the changes in serai-major axis, velocity, and 
period per orbit of a spacecraft under drag. In all of 
these equations, ballistic coefficient (B) was substituted 
for M/CjS r;r . The program was run using both radius and semi- 
major axis for the term "a" in the above equations. Results 
should be very similar since semi-major axis and radius are 
assumed to be equal in the formulation of the analytic 
equations and should be relatively equal for the small 
changes expected. Percent errors were calculated for the 
computed vs analytical results for changes in semi-major 
axis (SMA) , velocity (VEL) , and period (PER) , as well as for 
the percent differences between computed and analytical 
values of semi-major axis, velocity, and period after each 
orbit. Table I shows the results when radius is used in the 
analytical equations for "a", while Table II shows the 
results when semi-major axis is used in the analytical 
equations for "a". As shown in Table I, the percent error 
for the changes in semi-major axis ( %ASMA error) and period 
(%APER error) start small after one orbit, grow through the 
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Table I COMPUTED AND ANALYTIC ERRORS USING RADIUS 



ORBIT 

rr 


%ASMA 

ERROR 


%SMA 

ERROR 


%AP 

ERROR 


% PER 
ERROR 


%AV 

ERROR 


% VEL 
ERROR 


i 


0.4317 


. 00001 


0.4318 


. 00000 


7 .2875 


. 00015 


2 


1.6709 


. 00003 


1 . 6710 


. 00001 


5 . 1534 


. 00029 


3 


3 . 5508 


. 00001 


3 . 5509 


. 00002 


1.3233 


. 00043 


4 


5 .8015 


. 00005 


5.8017 


. 00005 


2 . 2979 


. 00056 


5 


8.0777 


. 00006 


8.0779 


. 00008 


6 . 6409 


. 00069 


6 


10.017 


. 00007 


10.017 


. 00012 


10.484 


. 00081 


7 


11.323 


. 00009 


11.326 


. 00015 


13 . 155 


. 00094 


8 


11 .860 


. 00010 


11.860 


. 00018 


14.265 


.00107 j 


9 


11.665 


. 00011 


11.666 


. 00020 


13.859 


| 

. 00120 


10 


10.948 


. 00012 


10.948 


. 00021 


12.377 


. 00134 



Table II COMPUTED AND ANALYTIC ERRORS USING SEMI-MAJOR AXIS 



ORBIT 

# 


%ASMA 

ERROR 


% SMA 
ERROR 


%AP 

ERROR 


% PER 
ERROR 


%AV 

ERROR 


% VEL 
ERROR 


1 


9.5531 


. 00001 


9.5531 


. 00002 


7 .2875 


.00015 ! 


2 


9.5531 


. 00003 


9.5532 


. 00005 


5.1534 


. 00029 


3 


9.5531 


.00004 


9.5532 


. 00006 


1.8283 


. 00043 


4 


9.5532 


. 00005 


9.5533 


. 00007 


2.2979 


. 00056 


5 


9.5532 


. 00006 


9.5534 


. 00009 


6.6409 


. 00069 


6 


9.5532 


. 00007 


9.5535 


. 00011 


10.484 


. 00081 


7 


9.5533 


. 0000® 


9.5535 


. 00013 


13 . 155 


. 00094 


8 


9 . 5533 


. 00010 


9.5536 


. 00015 


14 . 265 


.00107 


9 


9.5533 


. 00011 


9.5537 


. 00017 


13.859 


. 00120 


10 


9.5534 


. 00012 


9.5537 


. 00019 


12 . 377 


. 00134 
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fifth orbit, and then steady out at roughly eleven percent. 
This final error percentage closely corresponds with the 
relatively constant 9.5% error found when using semi-major 
axis in Table II. However, in both cases, as shown in 
Tables I and II, the percent error between the calculated 
and analytic semi-major axis (% SMA error) and period (% PER 
error) after each orbit is less than 0.0003%. The absolute 
percent error in the change in velocity (%AV error) and the 
difference between the calculated and analytic velocity (% 
VEL error) after each orbit are the same whether using semi- 
major axis or radius for "a", as shown in Tables I and II. 
Again, while the percent error in the change in velocity 
grows to greater than ten percent, the percent error between 
calculated and analytic velocity after each orbit very small 
(less than 0.002%). The reason for the ten percent errors 
between parameter changes appears to be the accuracy 
involved when the changes are so small (i.e. semi-major axis 
is changing by approximately 200 meters per orbit out of an 
initial orbit of 6638200 meters) . Thus an error of 10% in 
change of semi-major axis is only 20 meters, which is 
insignificant in terms of the initial radius. The analytic 
and computed results are thus in agreement and the drag 
portion of the program is validated. 
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3. Validation of thrust 

The thrust portion of the program was validated 
throughout the study by having drag present and the 
spacecraft maintain an orbital band by reboosting using 
thrusters . 

C. ORBIT CONTROL STRATEGY 

It is desired to design a control strategy that is 
simple and will maintain the spacecraft within a pre- 
specified radial band. Figure 2 shows a typical band, 
defining limits of radius and specific energy. 




Figure 2 Orbital band with radius and energy 
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All control strategies were run using numerous combinations 
of parameter values (i.e thrust vector angle, a, was varied 
from 0-90° in 5° increments, radial band was varied using 2, 
25 and 100 km bands, thrust was varied from 25 to 300 N., 
etc.)- For consistency, all cases presented here were run 
with the following baseline values 

R0= 6638 . 2km 
M- 2 000 0kg 

ThzusC = 3QQN 

Ballisticcoeffi cient -B = 15 Qkg/m z 
p 0 =9 . 4e-10kgr/'/7? 3 

(3 =2 . I2e-5km~ 1 
- sp = 3 OOsec 

a=70° 

Radialband= 2km 

which produce representative results corresponding to 

S=25000 

r=o . 060 

l. Control using radius 

Initial attempts were made to control the radial 
band using radius as a parameter, i.e. turn thrusters on 
when the radius drops to the low end of the band (RMIN) and 
turn them off when the satellite is reboosted to the high 
end of the band (RMAX) . As shown in Figures 1 of Appendix 
C, this control strategy was ineffective because the 
specific energy of the orbit continually increased and never 
leveled off (see Appendix C, Figure 2) . The eccentricity of 
this orbit grew to approximately 0.1. 



21 



2. Control using total energy 

Control was next attempted using total energy. The 
total energy of the initial orbit was calculated, as were 
the total energies of the spacecraft in circular orbits at 
the minimum and maximum of the radial band (see Figure 2 ) 
using 

(45) 

The thrusters were turned on when the spacecraft total 
energy dropped below the minimum orbit total energy and 
turned off when the spacecraft total energy rose above 
maximum orbit total energy. This method failed to control 
radius, as the spacecraft orbit always eventually decayed, 
but successfully controlled energy, as shown in Figures 3 
and 4 of Appendix C. While this method succeeded in 
maintaining a pre-determined radial band, the orbital radius 
was not maintained and decayed. 

3. control using radius and total energy 

A combination control strategy utilizing both radius 
and total energy was attempted. Thrusters were turned on 
when the orbital radius dropped to the minimum radius and 
turned off when total energy increased to the initial total 
energy. Figures 5 and 6 of Appendix C show the unsuccessful 
results. Total energy would seem to be the parameter of 
choice for maintaining the orbit vice specific energy. 
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However, if the total energy is maintained at a constant 
value the orbital radius must decrease and/or velocity must 
decrease since the spacecraft mass decreases due to 
thrusting. This can be shown using equation 13. 

4. Control using specific energy 

The next attempt at control involved spacecraft 
specific energy. This control was similar to the attempt 
using total energy, only specific energy was used. Specific 
energy is a function of orbital radius (for circular 
orbits) , so "conserving" specific energy should result in a 
"conservation" of radius (i.e. a constant average radius). 

As shown in Figure 7 of Appendix C, the orbital radius was 
maintained, however the orbital band was not. Conserving 
spacecraft specific energy, as shown in Figure 8 of Appendix 
C, failed to maintain the initial radius of the spacecraft. 

5. Control using radius and specific energy 

Control using radius and spacecraft specific energy 
is successful in maintaining radius and a radial band, as 
shown in Figures 9 and 10 of Appendix C. Thrusters are 
turned on based on the spacecraft orbit decaying below a 
specified radial minimum and spacecraft specific energy 
being below the initial specific energy. Thrusters remain 
on until specific energy increases above the initial value. 
This control strategy results in an orbit whose radius 
oscillates about the initial radius, but not necessarily at 
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the pre-specif ied value. The variation of the actual radial 



band from 
detail in 



the desired band varies and is examined in greater 
the next chapter. 



24 



V. 



ANALYSIS AND RESULTS 



Analysis of the validated control strategy is performed 
in five steps. A set of parameter values is chosen as the 
baseline case about which they will be varied. These values 
are selected as representative of a space-station/platform 
type satellite. The baseline values are 

i?=663 8 . 2km 
M= 20000kg 
J sp =300sec 
T=3 0 ON 

p 0 =9 . 4e-10icg r /.m 3 

B= 150kg/m 2 
1=25000 
P =2 . 12e-5 km' 1 

The simulation program is run for 100 orbits, with the 
equations of motion updated 5000 times per orbit, and output 
sampled 100 times per orbit. 

A. VARIATION OF THRUST ANGLE, a 

First, numerous cases are studied to determine the 
validity of an optimal thrust angle, a, of 70°. The program 
is run with combinations of variations of parameters about 
the baseline values, while varying a from 0-90° in 5° 
increments. In all cases, thrust angles of 65 and 70° are 
the only values that maintain a constant radial band. 

Figures 1-4 of Appendix D show results over the range of 60- 
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75° for the baseline parameters, and are representative of 
all cases. Values of a below 60° and above 75° produce 
results that degrade significantly from those shown. Semi- 
major axis (SMA) and eccentricity are also compared 
throughout the analysis for determining optimal values of a. 
Figures 5 and 6 of Appendix D show results consistent with a 
maintained radial band (a=70°) , while Figures 7 and 3 show 
results consistent with a radial band not being maintained 
(a=50°) . As can be seen, eccentricity must level off to a 
constant value to maintain a radial band, while SMA alone is 
not an indicator of orbital maintenance. 

The mass of fuel required to maintain orbit generally 
increases as a increases from 0 to 90°, with some minor 
local fluctuations. Table III shows the mass of fuel burned 
in kilograms (kg) to maintain 2, 25, and 100 km pre- 
specified bands, with a varying from 60-75° and baseline 

Table III MASS OF FUEL BURNED (kg) 



Band (km) 


a=60° 


a=65° 


a=7 0° 


a=75° 


2 


2643 


2013 


2174 


2683 


25 


2643 


2828 


3243 


3595 


100 


3020 


2853 


2826 


3787 



values for other parameters, and is typical of the trend 
seen. Thus, from a fuel usage standpoint a should be 
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minimized, but from a radial band standpoint a must be in 
the range of 65-70°. 

This shows that a constant thrust angle, a, of 
approximately 65-70° is necessary to maintain a radial band. 
The following analyses are conducted with a range of 60-75° 
for a to minimize the data collection while ensuring the 
continuing validity of this result. 

B. VARIATION OF PRE-SPECIFIED RADIAL BANDWIDTH 

Once a is chosen, it is desired to know how the pre- 
specified radial bandwidth affects the results. While 
maintaining a within the 60-75° range and using the baseline 
values for other parameters, the program is run using 
specified radial bandwidths of 1,2,5,10,25,100, and 200 km. 
Although the large bandwidths are not practical, they are 
used to determine the existance of any limiting conditions. 
These results, shown in Table IV (first entry is fuel burned 
in kg, second entry is actual radial band maintained in km) , 
show that fuel required to maintain the orbit generally 
increases as the controlled radial bandwidth increases, with 
the exception being the 100 km band. Although the error 
between the radial band maintained and the pre-specif ied 
band is minimized as a increases, the maintained band cannot 
be correlated to the pre-specif ied band (changing any 
parameter changes the maintained band) . For this reason, 
the radial bands discussed throughout this analysis are pre- 
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Table IV MASS OF FUEL BURNED (kg) AND ACTUAL RADIAL BAND 
MAINTAINED (km) 



Band 


a=60 


a=65 


a=7 0 


a=75 


! i 


1537/0.7 


1684/0.7 


2025/1.0 


2627/1.2 


2 


2049/2 


2013/1.8 


2174/1.8 


2683/2 


5 


2516/12 


2628/8 


2583/6 


3065/5 


10 


2392/25 


2554/16 


3191/13 


3602/11 


25 


2643/55 


2828/43 


3243/35 


3595/30 


50 


2812/90 


3049/75 


3136/65 


4289/58 


100 


3020/130 


2859/125 


2826/60 


3787/75 


200 


2993/130 


3549/140 


4470/120 


5743/135 



specified bands. Figures 1-3 of Appendix E show the mass of 
fuel burned using bang-bang control and forced Keplerian 
control for 2, 25, and 100 km bands with a=70°. As can be 
seen, the fuel used increases in a generally linear manner 
as the number of orbits increases. For the 100 km band, the 
fuel used is linear when the points at the end of each burn 
are connected, as indicated by the dotted line in Figure 3. 
Thus, these results can be extrapolated for any number of 
orbits to give fuel usage vs Keplerian usage. However, the 
thrusters do not fire at the same point in each orbit (the 
thrusters fire approximately every three-fourths of an orbit 
for the cases of a small radial band and, thus firing 
precesses each orbit) . These results are consistent 
throughout the analysis. 
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A constant radial bana is maintained for a of 



approximately 65-70° but not for 60 and 75°, as shown in 
Figures 4-7 of Appendix E, which again verifies the required 
thrust angle. These figures are for a pre-specif iea raaial 
bandwidth of 25 km, the results are similar for banawidths 
of 1-50 km. Figure 8 of Appendix E shows the results for a 
pre-specif ied bandwidth of 100 km and an a of 70° (this is 
representative of results for a=60-75° and bands of 100 km 
and greater). As can be seen, the band is wide enough that 
thrusters fire infrequently to reboost the satellite. The 
actual radial bandwidth maintained was then compared to the 
specified band. Based on these results, all subsequent 
analysis is performed for three cases of radial bandwidth; 

2, 25, and 100 km. This will give results for a small, 
medium and large bandwidths, while examining the unusual 
results at 100 km. 

C. VARIATION OF ATMOSPHERIC SCALE HEIGHT, 0 

Of the three variable parameters in the non- 
dimensionalized equations, /3 should have the most 
predictable effect on the results and is analyzed first. 
Using baseline values for other parameters, /3 is varied from 
1.0e-10 to l.0e-2, and results are examined for logical 
trends. Logic would dictate that increasing the value of 0 
will increase the rate of orbital decay within the band, 
requiring more thruster firings to reboost the satellite, 
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and thus require more fuel. This is indeed the case as the 
results in Appendix F show. Figure 1 shows the fuel used as 
a function of j3 for a pre-specif ied 2 km band, while Figures 
2 and 3 show the same results for 25 and 100 km bands, 
respectively. All three plots are similar in appearance 
indicating changes in affect orbital decay uniformly over 
the range of bands. Below a certain value, le-5 for 2 km 
band, changes in (3 have no effect on fuel usage. Above a 
certain value, le-3 for 2 km band, the orbit decays too 
quickly for thrusters to control. In between these values, 
increases in increasingly affect the fuel used. Since 
changes in (3 do not affect the forced Kepler ian solution, 
the mass of fuel required for forced Keplerian maintenance 
for all cases is approximately 649 kg. Thus, smaller values 
of /3 allow the bang-bang control to be closer in efficiency 
to forced Keplerian motion. Analysis of the plots for 
radial band again show that a radial band is not maintained 
for 60 or 75° thrust angle, but is maintained and optimized 
in the 65-7 0° range. From these results, it is shown that 
changes in (3 affect the efficiency of bang-bang control as 
expected and are predictable. 

D. VARIATION OF THRUST, T 

The next variable to be analyzed is thrust (T) . As 
explained in Chapter 3, either thrust or nondimensionalized 
thrust can be varied and the results will be the same. 
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Thrust is chosen because variations in thrust are more 
easily related to and understood. For this analysis, thrust 
is varied over the range 10-1000 N, while all otner 
parameters are fixed at the baseline value. Figures 1-3 of 
Appendix G show the fuel required, for a=70°, to maintain 2, 
25, and 100 km pre-specif ied bands, respectively. 

Trends are similar on all three graphs. Fuel usage 
starts off low for small values of thrust, with thrust being 
small enough to produce quasi-forced Keplerian motion 
(thrusters fire for majority of time to reboost satellite) . 
Forced Keplerian thrust is determined to be approximately 
3.3 N for the reference conditions. To obtain this value of 
transverse thrust using a thrust angle of 70° would require 
a thrust of approximately 11.1 N. Thus, small values of 
thrust may produce a quasi-forced Keplerian motion with 
respect to the transverse thrust component. This is 
followed by a period of increased fuel usage as thrust 
increases. In the 2 km band case, increasing thrust further 
shows a drop in fuel usage followed by a relatively constant 
fuel usage as thrust increases above 300 N. The 25 km band 
case shows more fluctuation in fuel usage as thrust 
increases, followed again by a relatively constant fuel 
usage as thrust increases above 300 N. For the 100 km band 
case, increasing thruster size generally increases fuel 
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usage except for a fairly constant region in the range of 
50-250 N. 

Thus, other than a possible minimizing of fuel usage due 
to this apparent quasi-forced Keplerian motion, changes in 
thrust appear to have little effect on the optimality 
problem. The smaller fuel usage totals at low thrust (25-50 
N) may be indicative of approaching a quasi-forced Keplerian 
motion (except thruster angle is 70° vice 0°) . For all three 
cases, as thrust increases, fuel usage appears to approach a 
constant value. An optimal thrust angle of 65-70° is again 
verified in this portion of the analysis. 

E. VARIATIONS IN NONDIMENS I ONAL I ZED BALLISTIC COEFFICIENT 

The last parameter to be analyzed is nondimensional 
ballistic coefficient, "B. By nondimensionalizing, 
variations in B will encompass variations in the individual 
parameters of mass, atmospheric density, coefficient of 
drag, and reference surface area, as well as combinations of 
variations involving these parameters. 

Analysis is conducted by varying B over the range 4e3- 
2e5 and comparing mass of fuel burned using bang-bang 
control (MFB) to mass of fuel burned using forced Keplerian 
motion (MFBK) . This analysis is conducted for the three 
bands of 2, 25, and 100 km using baseline values for other 
parameters while varying a from 60-75°. Figures 1 and 2 of 
Appendix H show the results for the 2 km band case, results 
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for rhe 25 and 100 km bands are similar. As seen, the mass 
of fuel burned for both bang-bang control and forced 
Kepierian motion decreases as B increases, approaching a 
constant value as B increases above 2e5. For optimality, 
maximizing B minimizes fuel usage. 

However, this does not give a feel for how efficient 
bang-bang control is compared to forced Kepierian motion. 

To analyze this, the ratio of MFB to MFBK is plotted versus 
B for 2, 25, and 100 km bands and a of 65 and 70°. These 
results are shown in Figures 3-8 of Appendix H. Again, B is 
varied over the range 4e3-2e5. Values of B below and above 
this range (and the range plotted for 100 km band) cause the 
orbit to decay too rapidly or too slowly, respectively. 
Logically, this is analogous to increasing the atmospheric 
density to a value that causes the orbit to decay faster 
than thrusters can maintain, or decreasing the atmospheric 
density to a value that causes the orbit to decay so slowly 
that the thrusters never fire in 100 orbits (thus giving no 
data since at least one firing of the thrusters in bang-bang 
control is required for comparison to forced Kepierian 
motion) . 

As seen in Figures 3-8, the general results for a of 65 
and 70° are similar for all three bands. Also, the results 
for all three bands show a roughly constant value for the 
ratio of MFB to MFBK, with the 2 km band being most 
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efficient compared to forced Keplerian motion, followed by 
the 100 km band and then the 25 km band. A value of 1 for 
the ratio would indicate chat the bang-bang control was 
equally efficient to Keplerian motion. Values less than 1 
would indicate that the bang-bang control is superior to 
Keplerian motion and would numerically validate the 
theoretical results of Ross and Melton [Ref. 3]. The 
results shown are consistently above a ratio of 3, 
indicating that forced Keplerian motion is significantly 
more efficient than bang-bang control over all values of B. 

Throughout this analysis, the optimality of a in the 
range of 65-70° was verified. This result has remained 
valid throughout all of the analyses conducted, indicating 
that this is a globally valid result. 



34 



VI. CONCLUSIONS AND RECOMMENDATIONS 



The purpose of this thesis was to determine if there are 
any cases where fixed angle, bang-bang control, orbital 
maintenance approaches or exceeds the efficiency of forced 
Keplerian motion as proposed by Ross and Melton [Ref. 3]. 

By nondimensionaliz ing the equations, a thorough study of 
the effects of all parameter combinations was possible. 
Analysis of the results indicate that, in all cases, forced 
Keplerian motion is superior to bang-bang control. There 
are no cases where bang-bang control even remotely 
approaches the efficiency of forced Keplerian motion. Thus, 
these results appear to be globally valid. 

The optimality question for thrust angle also appears to 
be resolved globally. In all cases, the thrust angle is 
required to be between 65 and 70° in order to maintain a 
constant radial band. Generally, minimizing this angle 
improves fuel consumption, so the angle should be minimized 
as far as practicable while ensuring a constant radial band. 

Maintenance of a pre-specif ied radial bandwidth is not 
always possible or predictable for this control strategy. 

The error between specified and maintained band can be 
minimized through proper bandwidth selection. Selecting a 
relatively small, 2 km, or large, 100 km, band minimizes 



35 



the error. Also, maximizing the thrust angle, a, generally 
minimizes the error, and thus, a should be set as high as 
possible, while maintaining a constant radial band. Since 
this contradicts the requirement to minimize a for fuel 
savings, trade-offs must be considered. Minimizing the 
specified bandwidth moves the efficiency of bang-bang 
control closer to that of forced Kepler ian motion. 

The effects of changes in atmospheric scale height, ( 3 , 
are logical and predictable on the effects of the mass of 
fuel burned in bang-bang control. Since, by definition, 
changes in (3 do not affect forced Kepierian results, the 
changes in mass burned in bang-bang control directly 
correlate to changes in efficiency. Thus, the lower the 
value of /?, the more efficient bang-bang control is in 
relation to Kepierian motion. 

The effects of varying thrust are not so predictable; 
however, general trends can be inferred from the results. 

Low values of thrust appear to approach quasi-forced 
Kepierian motion and are more efficient. Higher values of 
thrust increase the mass of fuel burned, however, a maximum 
limit appears to be approached. As in the case of ( 3 , 
changes in thrust do not affect the mass of fuel burned in 
Kepierian control, thus the mass of fuel burned in bang-bang 
control is indicative of the efficiency. 

This is not the case for changes in B, where changes 
affect the final Kepierian fuel usage. Thus, in order to 
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examine efficiency, the ratio of mass of fuel burned using 
bang-bang control to mass of fuel burned with Keplerian 
motion must be used. It was initially expected that cnanges 
_n B would have noticeable effects on efficiency and would 
show a trend that would indicate values where bang-bang 
control efficiency approached Keplerian efficiency. This is 
not the case, as changes in B seem to have little effect on 
the ratio of efficiencies. 

Thus, the superiority of forced Keplerian motion to 
fixed angle bang-bang control is globally confirmed. 

However, further research should be conducted into the 
solutions proposed by Ross and Melton [Ref. 3]. Control 
using variable thrust angle bang-bang control (so-called 
Primer-vectoring) should be examined for comparison to these 
results. Also, fixed angle bang-bang control, as described 
here, should be compared with other innovative methods of 
control, as well as with the "Lambert control" used with the 
shuttle, for possible fuel savings for space-station type 
platforms . 

The results of this thesis also suggest areas for 
further study. A method for maintaining a pre-specif ied 
band is highly desirable, but was not achieved in this 
study. Although the global validity of the thrust angle 
being required to be between 65 and 70° was shown, the 
reason for this has not been studied. The problem of 
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micro-gravity constraints and its .effect on thruster firings 
and control strategies has also nor been studied, but would 
be a critical factor in any control design. The facr that 
forced Keplerian motion has been shown to be on the order of 
three times more efficient than bang-bang control indicates 
that there is much room for research into further 
optimization techniques for bang-bang control. Finally, the 
question of practicality has not been addressed in relation 
to both forced Keplerian motion (generally not very 
practicable but very efficient) and bang-bang control (very 
practicable but not very efficient) . 
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APPENDIX A 



A. PROGRAM LISTING 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



OBJECTIVE: COMPARISON OF FORCED KEPLARIAN TO BANG-BANG 
CONTROL OF ORBITAL ALTITUDE 



VARIABLE DEFINITIONS 

XBAR ( 1 ) =NONDIMENSIONALIZED ORBITAL RADIUS 

XBAR ( 2 ) =NONDIMENSIONALI ZED ORBITAL RADIAL VELOCITY 

XBAR ( 3 ) =THETA 

XBAR ( 4 ) =NONDIMENSIONALI ZED ORBITAL ANGULAR VELOCITY 
XBDOT ( 1 ) =DERIVATIVE OF XBAR ( 1 ) 

XBDOT ( 2 ) =DERI VATIVE OF XBAR ( 2 ) 

XBDOT ( 3 ) =DERIVATIVE OF XBAR ( 3 ) 

XBDOT( 4 )=DERI VATIVE OF XBAR ( 4 ) 

R0= INITIAL ORBITAL RADIUS 
R=ORBITAL RADIUS 

RREF=REFERENCE ALTITUDE FOR DENSITY 
D=DRAG (N) 

DK=DRAG KEPLERIAN 
E=SPECIFIC ENERGY 

SO=SPECI FIC ENERGY AT INITIAL RADIUS RO 
M=S/C MASS (KG) 

MO=INITIAL S/C MASS (KG) 

ME=ARB I TRARY MASS FOR NONDIMENSIONALIZING (20000KG) 

MBAR=NONDIMENSIONALI ZED MASS 

MF=MASS OF FUEL BURNED IN TIME INCREMENT 

MFK=MASS OF FUEL BURNED PER TIME INCREMENT KEPLERIAN 

MFT=MASS OF FUEL BURNED TOTAL 

MFTK=MASS OF FUEL BURNED TOTAL WITH KEPLERIAN MOTION 
GAMMAR= FLIGHT PATH ANGLE (RAD) 

GAMMAD= FLIGHT PATH ANGLE (DEG) 

TH=THRUST (N) 

THK=KEPLERIAN THRUST 

THMAX= BLOWDOWN (MAX) THRUST (N) 

THBAR=NONDIMENSIONALI ZED THRUST = ( TH*TAU* * 2 ) / ( MR ) 
THBM=NOND IMENSIONALI ZED BLOWDOWN (MAX) THRUST 
ALPHAR=THRUST ANGLE (RAD) 

B=BALLISTIC COEFFICIENT (M/(CD*S)) 

BBAR=NOND IMENSIONALI ZED B =B/(RHOO*RE) 

RHO=CALCULATED ATMOSPHERIC DENSITY 
RHOO=REFERENCE ATMOSPHERIC DENSITY 
SPI= SPECIFIC IMPULSE 
V=VELOCITY 

VBAR=NOND IMENSIONALI ZED VELOCITY 
G=G RAVITATIONAL ACCELERATION 
SMA=SEMI -MAJOR AXIS 
ECC=ECCENTRICITY 
T=TIME 

T I NC= INCREMENT OF TIME (STEP SIZE) 
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non noi-*nno noooo ooo non ooo 



TF=END TIME 

PTI=PRINT TIME INTERVAL 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



CONSTANTS 

RE=RADIUS OF EARTH 

MU=EARTHS GRAVITAIONAL CONSTANT 

TAU=ORB I TAL PERIOD AT RE 



START PROGRAM 
PROGRAM ORBMAINT 



VARIABLE DECLARATION 
IMPLICIT REAL* 8 (A-H,M-Z ) 
DIMENSION XBAR ( 4 ) , XBDOT ( 4 ) 



CONSTANT DEFINITIONS 
PI =DATAN ( 1 . 0D+00 ) *4.0D+00 
MU=3 .98601208133D+14 
G=9 . 806D+0 
RE=6 . 3782D+6 

TAU=2 . *PI *RE* * 1 . 5 /MU* *0 . 5 
J=4 

ME=20000 . D+0 

INITIALIZE FUEL USAGE TOTALS 

MFT=0. OD+O 
MFTK=0 . 0D+0 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

MAIN PROGRAM 



OPEN ( 10 , FILE= ' inp ' , STATUS= ' OLD ' ) 

OPEN ( 20 , FILE= ' outp ' , STATUS= ' NEW ' ) 

OPEN ( 30 , FILE= ' orbp ' , STATUS='NEW' ) 

OPEN ( 40 , FILE= ' ratp ' , STATUS= ' NEW' ) 

READ INPUTS AND INITIALIZE PARAMETERS 

READ ( 10 , 1 ) R0 , M0 , THMAX , TB , TFB , TINCB , PT IB , BBAR , BETA , SP I 
FORMAT (10(/,21X,D13.7) ) 



M=M0 

MBAR=M/ME 

TAUR0=2. *PI*R0**1.5/MU**0.5 
TINC=TINCB*TAU 
THBM=THMAX*TAU*TAU/ (ME*RE) 

READ PARAMETERS THRUST ANGLE, DESIRED BAND, THRUST FACTOR 
c PRINT*, 'ENTER ALPHA' 



40 



c 

c 



c 

c 

c 



c 

c 

c 



c 

c 

c 

c 



c 

c 

c 



c 

c 

c 



c 



c 

c 

c 



c 



c 

c 



READ* , ALPHA 

ALPHAR=ALPHA*PI / 180 . 0D+0 



PRINT* , ' ENTER BAND , { KM ) ' 

READ* , BAND 
3AND=BAND * 1000 . 0D + 0 



PRINT* , 7 ENTER THFAC ' 

READ* , THFAC 
THFAC=THFAC* 1 . 0D+0 

INITIALIZE COUNTERS USED IN PROGRAM 

INDEX=0 
KOUNT= 1 
TRAD=0 . 0D+0 



INITIALIZE NONDIMENSIONALIZED STATE VARIABLES 

XBAR ( 1 ) =R0 / RE 
XBAR ( 2 ) =0D+0 
XBAR ( 3 ) = 0D + 0 

XBAR ( 4 ) =2 . 0D+0*PI *TAU/TAUR0 
COMPUTE INITIAL ORBITS VELOCITY AND ENERGY 



V0= ( MU/RO ) * *0 . 5 
E0=( (V0*V0) /2. -MU/RO) 

COMPUTE RADIUS , VELOCITY , ENERGY FOR DESIRED BAND MIN AND MAX 



RMIN=RO-BAND/2.0 

VMIN=(MU/RMIN) **0.5 

EMIN= ( ( VMIN*VMIN ) /2 . -MU/RMIN ) 

RMAX=R0+BAND/2.0 

VMAX= ( MU/RMAX ) * *0 . 5 

EMAX=( ( VMAX*VMAX ) / 2 . 0-MU/RMAX ) 



SET UP OUTPUT FILE HEADERS 



IBAND=INT ( BAND ) 
ITHMAX= INT ( THMAX ) 
IALPHA= INT (ALPHA) 



WRITE ( 20 , * ) ' ! 
WRITE (20,*) ' ! 
WRITE ( 30 , * ) ' ! 
WRITE (40 , * ) ' \ 



! BAND= ' , I BAND , ' 
! EMIN= ' , EMIN, ' 

! BAND= ' , I BAND, ' 
i BAND= ' , I BAND , ' 



THMAX= ' , ITHMAX, ' 
EMAX= ' , EMAX 
THMAX= 9 , ITHMAX, # 
THMAX= # , ITHMAX, # 



ALPHA= # , I ALPHA 

ALP HA = ' , I ALPHA 
ALPHA= ' , I ALPHA 



WRITE (20, * ) ' ! 1 
ANGM 


ORBITS RADIUS 

ENERGY ' 


VELOCITY 


MFB 


MFBK 


WRITE (20,*) ' ! ! 


(KM) 


(KM/SEC) 


(KG) 


(KG) ' 


WRITE ( 30 ,*)'!! 
PERIOD ' 


ORBITS SMA 


ECC 


APOGEE 


PERIGEE 


WRITE (40, * ) ' ! ! 

GAMMA' 


TIME ORBITS 


DRAG 


THRUST 


MASS 
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c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

100 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 



c 



CALCULATIONS 

CALL DRAG TO SET INITIAL VALUE FOR FORCED KEPLERIAN MOTION 
USED TO COMPUTE KEPLERIAN THRUST 

CALL DRAG ( BBAR , RE , XBAR , BETA , MBAR , VBAR , D ) 

DK=D 

TK=DK*MBAR*ME *RE / ( TAU*TAU ) 

CALL SUBROUTINE DRAG TO COMPUTE PRESENT DRAG 

CALL DRAG ( BBAR , RE , XBAR , BETA , MBAR , VBAR , D ) 

CALCULATE E TO SEND TO THRUST SUBROUTINE TO DETERMINE 
IF THRUSTER SHOULD BE ON OR OFF 

E=( RE *VBAR/TAU) **2.0/2. OD+O-MU/ (XBAR(l) *RE ) 

CALL SUBROUTINE THRUST TO DETERMINE IF THRUSTERS ARE ON OR 
OFF AND THE SET THE VALUE OF THE THRUST TERM 

CALL THRUST ( EO , XBAR , THBAR , THFAC , RMIN , EMAX , E , ETH , THBM ) 

CALL SUBROUTINE EQN TO UPDATE ORBITAL EQUATIONS 

CALL EQN ( D , XBAR , XBDOT , VBAR , THBAR , ALPHAR , MBAR , PI ) 

CALL SUBROUTINE RK4 TO USE RUNGE-KUTTA FOR INTEGRATION 
OVER TIME 

CALL RK4 ( TB , XBAR, XBDOT , J, TINCB , INDEX) 

CHECK TO SEE IF RUNGE-KUTTA COMPLETE 
IF (INDEX .NE. 0) GO TO 100 



RUNGE-KUTTA DONE FOR THIS ITTERATION OF ORBIT , 
UPDATE VARIABLES FOR CALCULATION OF ORBITAL 
PARAMETERS AND FOR OUTPUT 

TH=THBAR*ME*RE/ ( TAU*TAU ) 

R=XBAR ( 1 ) *RE 
V=RE*VBAR/TAU 
MF=TH*TINC/ (SPI*G) 

M=M-MF 

MFT=MFT+MF 

GAMMAR=ATAN ( XBAR ( 2 ) / ( XBAR ( 1 ) *XBAR ( 4 ) ) ) 
GAMMAD=GAMMAR* 180 . OD+O/PI 
ANGM=R*V*COS (GAMMAR) 

ENERGY=( (V*V) /2.0)-(MU/R) 

TRAD =TRAD +T I NCB 
OANG=TRAD*2 . *PI 

COMPUTE KEPLERIAN FUEL VALUES 

MFK=TK*TINC / ( SPI *G ) 

MFTK=MFTK+MFK 
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C CALL SUBROUTINE ORBPAR TO COMPUTE ORBITAL PARAMETERS 

C 

CALL ORBPAR ( ENERGY , ANGM , R , MU , PI , ECC , SMA , APOGEE , PERIGE , PERIOD ) 

C 

C COMPUTE NEW MBAR AND PUT R,V IN OUTPUT FORM 

C 

MBAR=M / ME 
R=R/ 1000 . 0 
V=V/1000.0 

c 

c 

C DATA OUT 

C 

C CHECK TO SEE IF TIME TO PRINT OUTPUT 

C 

IF ( KOUNT .LT. DNINT ( PTIB/TINCB ) ) GO TO 200 
C 

C PRINT OUTPUTS TO OUTPUT FILES OUT , RAT , ORB 

C 

WRITE ( 2 0 , 2 ) TB , R , V , MFT , MFTK , ANGM , E 

2 FORMAT (2X,F6.2,2X,F9.3,2X,F7.4,2X,F8.2,2X,F8.2,2X,F12.0, 

*2X , F10 . 0 ) 

C 

WRITE (30,3) TB , SMA , ECC , APOGEE , PERIGE , PERIOD 

3 FORMAT (2X,F7. 2, IX, F10 . 3 , 3X , F4 . 3 , 3X, F10 . 3 , IX , F10 . 3 , 2X , F8 . 4 ) 

C 

WRITE ( 40 , 4 ) T , TB , D , TH , M , GAMMAD 

4 FORMAT (2X,F8.0,1X,F6.2,2X,F12.9,3X,F5.0,5X,F9.3, IX, F10.4) 

C 

C RESET COUNTER 

C 

KOUNT=0 

C 

C CHECK TO SEE IF ORBIT COMPLETE, RESET ORBITAL ANGLE 

C 

IF (TRAD .GE. 1.0D+0) TRAD=0.0D+0 
C 

C UPDATE COUNTER 

200 KOUNT=KOUNT+1 
C 

C CHECK FOR END OF PROGRAM 

C 

IF ( TB .LT. TFB ) GO TO 100 
C 

END 

C 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c THIS SUBROUTINE CALCULATES DRAG FOR USE IN THE EQUATIONS 

C OF MOTION. 

C 

SUBROUTINE DRAG ( BBAR , RE , XBAR , BETA , MBAR , VBAR , D ) 

IMPLICIT REAL*8(A-H,M-Z) 

C 

DIMENSION XBAR ( 4 ) 

C 

C SET REFERENCE ALTITUDE FOR ATMOSPHERIC DENSITY CALCULATIONS 

C 

RREF=6 . 63814 5D+6 
C 
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c COMPUTE PRESENT NONDIMENSIONALIZED VELOCITY 

VBAR=( ( XBAR ( 2 ) *XBAR ( 2 ) ) + ( XBAR ( 1 ) *XBAR ( 4 ) )**2)**0.5 
C 

C COMPUTE EXPONENTIAL TERM IN DRAG EQN 

C 

E=EXP ( -BETA*RE* ( XBAR ( 1 ) - ( RREF/RE ) ) ) 

C 

C COMPUTE DRAG TERM 

C 

D=E*VBAR*VBAR/ ( 2 . 0 *MBAR*BBAR ) 

C 

C 

RETURN 

END 

C 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c THIS SUBROUTINE UPDATES THE EQUATIONS 

C OF MOTION 

C 

SUBROUTINE EQN ( D , XBAR , XBDOT, VBAR , THBAR , ALPHAR , MBAR , P I ) 
IMPLICIT REAL *8 ( A-H , M-Z ) 

C 

DIMENSION XBAR ( 4 ) , XBDOT ( 4 ) 

C 

C 

A=XBAR ( 1 ) *XBAR ( 4 ) *XBAR ( 4 ) 

B=4 .0*PI*PI/ ( XBAR ( 1 ) *XBAR( 1 ) ) 

C=D* (XBAR(2) /VBAR) 

E=THBAR*SIN (ALPHAR) /MBAR 
C 
C 

F=2 . 0*XBAR ( 4 ) *XBAR ( 2 ) /XBAR ( 1 ) 

G= ( D / XBAR ( 1 ) ) * ( XBAR ( 1 ) *XBAR ( 4 ) /VBAR ) 

H=THBAR*COS (ALPHAR)/ (MBAR* XBAR(l) ) 

C 

C 

XBDOT ( 1 ) =XBAR ( 2 ) 

XBDOT ( 2 ) =A-B-C+E 
XBDOT ( 3 ) =XBAR ( 4 ) 

XBDOT ( 4 ) =-F-G+H 
C 
C 

RETURN 

END 

C 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C THIS SUBROUTINE CALCULATES ORBITAL PARAMETERS 

C 

SUBROUTINE ORB PAR (ENERGY , ANGM , R, MU , PI , ECC , SMA, APOGEE , PERIGE , 
*PERIOD ) 

IMPLICIT REAL* 8 (A-H, M-Z) 

C 

C PROBLEM COULD ARISE IN COMPUTING ECCENTRICITY IF VERY 

C SMALL- COULD TRY TO TAKE SQRT OF NEG NUMBER. TO 

C PREVENT, IF ECCENTRICITY IS LESS THAN IE-6, SET 

C EQUAL TO 0. 

C 
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?ROB= ( 1 . t-2 . * ENERGY *ANGM*ANGM/ ( MU*MU ) ) 

IF ( DABS ( PROB ) .LT. l.D-12) THEN 
ECC=0 . 0D+0 
ELSE 

ECC=PROB* *0.5 

ENDIF 

C 

SMA= ( -MU/ ( 2 . ^ENERGY ) ) / 1000 . 0 

C 

C 

SMAM=SMA* 1000.0 
APOGEE=SMA* ( 1 . 0+ECC ) 

PERIGE=SMA* ( 1.0-ECC) 

PERIOD =2 . 0*PI * (SMAM**1. 5 ) / (MU**0. 5 ) /60 .0 

C 

RETURN 

END 

C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C THIS SUBROUTINE DOES A FOURTH LEVEL RUNGE-KUTTA 

C INTEGRATION. 

C 

SUBROUTINE RK4 ( TB , XBAR , XBDOT , J , TINCB , INDEX ) 

IMPLICIT REAL * 8 ( A-H , M- Z ) 

INTEGER INDEX, I 
C 

DIMENSION XBAR ( 4 ) , XBDOT ( 4 ) , SAVED ( 4 ) , SAVEX ( 4 ) 

C 

INDEX=INDEX+1 
GO TO (1,2, 3, 4) , INDEX 
C 

1 DO 10 1=1, J 
SAVEX ( I ) =XBAR ( I ) 

SAVED ( I ) =XBDOT ( I ) 

10 XBAR ( I ) =SAVEX ( I ) +0 . 5D+0*TINCB*XBDOT ( I ) 

TB=TB + 0 . 5D+0 *TINCB 
RETURN 
C 
C 

2 DO 20 1=1, J 

SAVED ( I ) =SAVED ( I ) +2 . D+0*XBDOT ( I ) 

20 XBAR( I ) = SAVEX ( I ) +0 . 5D+0*TINCB*XBDOT ( I ) 

RETURN 

C 

C 

3 DO 30 1=1, J 

SAVED ( I ) = SAVED ( I ) +2 . D+0*XBDOT ( I ) 

30 XBAR ( I ) =SAVEX ( I ) +TINCB*XBDOT ( I ) 

TB=TB+0 . 5D+0*TINCB 
RETURN 
C 
C 

4 DO 40 1=1, J 

40 XBAR ( I ) =SAVEX ( I ) +TINCB / 6 . D+0 * ( SAVED ( I ) +XBDOT ( I ) ) 

INDEX=0 

RETURN 

C 

C 
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no Qoi-oon non non non non ono 



END 



C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C THIS SUBROUTINE DETERMINES IF THRUSTERS SHOULD BE ON OR 

C OFF TO TRY TO STAY WITHIN THE DESIRED BAND 

C 

SUBROUTINE THRUST ( EO , XBAR , THBAR , THFAC , RMIN , EMAX , E , ETH , 
*THBM) 

IMPLICIT REAL *8 ( A-H, M-Z ) 

C 

DIMENSION XBAR ( 4 ) , XBDOT ( 4 ) 

C 

RE=6 . 3782D+6 

C 

C COMPUTE NONDIMENSIONALIZED POINT FOR THRUSTERS TO TURN ON 

C BASED ON RADIUS. 

C 

RTH=RMIN/RE 

C 

c ETH=EMAX- ( EMAX-EMIN ) ^ THFAC 

C 

C IF THRUSTERS ARE ON GO TO SEE IF TIME TO TURN OFF 

C 

IF (THBAR .EQ. THBM ) GO TO 100 

CHECK IF RADIUS BELOW BAND 
IF ( XBAR ( 1 ) . LE . RTH ) THEN 

CHECK IF ENERGY BELOW INITIAL ENERGY 
IF (E .LE. EO) THEN 

TURN THRUSTERS ON 
THBAR=THBM 

ELSE 

OTHERWISE LEAVE OFF 
THBAR=0D+0 

END IF 

RADIUS NOT BELOW BAND, LEAVE OFF THRUSTERS 

ELSE 

THBAR=0D+0 

ENDIF 

IF THRUSTERS ARE ON, SEE IF TIME TO TURN OFF 
IF (E .GE. EO) THBAR=0D+0 



RETURN 

END 



CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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B. SAMPLE INPUT FILE 



S. 6381450D+06 
2 . OOOOOOOD+04 
3 . 0000000D+02 
0. OOOOOOOD+OO 
1. 0000000D+02 
2 . 0000000D-O4 
1. 0000000D-01 
2 . 5000000D+04 
2 . 12 00000D-05 
3 . OOOOOOOD+02 



C. SAMPLE OUTPUT FILE 



BAND = 


25000 


THMAX = 


300 


EMIN = 


-30080181. 4S670808 


EMAX = 


ORBITS 


RADIUS 


VELOCITY 


MFB 




(KM) 


(KM/SEC) 


(KG) 


1.00 


6636.385 


7.7501 


0.00 


2.00 


6634.580 


7.7512 


0.00 


3.00 


6632.737 


7.7523 


0.00 


4.00 


6630.855 


7.7534 


0.00 


5.00 


6628.930 


7.7545 


0.00 


6.00 


6626.951 


7.7557 


0.00 


7.00 


6634.414 


7.7513 


110.11 


8.00 


6628.660 


7.7600 


187.56 


9.00 


6625.232 


7.7641 


224.78 


10.00 


6622.517 


7.7673 


259.52 


11.00 


6620.638 


7.7695 


293.54 


12.00 


6619.667 


7.7706 


327.04 


13.00 


6619.631 


7.7706 


359.91 


14.00 


6620.517 


7.7696 


392.48 


15.00 


6622.261 


7.7676 


425.88 


16.00 


6624.814 


7.7642 


452.25 


17.00 


6628.234 


7.7599 


483.47 


18.00 


6632.050 


7.7555 


519.87 


19.00 


6636.064 


7.7509 


556.16 


20.00 


6640.091 


7.7462 


592.24 


21.00 


6643.950 


7.7417 


628.33 


22.00 


6647.468 


7.7377 


663.58 


23.00 


6650.501 


7.7342 


699.25 


24.00 


6652.923 


7.7314 


734.93 



ALPHA = 70 

-29967109.06483507 

MFBK ANGM ENERGY 

(KG) 

6.49 51432570538 . - 30031067 . 
12.98 51425887643 . - 30038873 . 
19.47 51418951743 . - 30046977 . 
25.96 51411742642 . - 30055404 . 
32.45 51404237633 . - 30064181 . 
38.94 51396411113 . - 30073338 . 
45.43 51425131788 . - 30039612 . 
51.91 51438552722 . - 30023833 . 
58.40 51438811678 . - 30023532 . 
64.89 51438814899 . - 30023533 . 
71.38 51438812800 . - 30023539 . 
77.87 51438816720 . - 30023537 . 
84.36 51438831447 . - 30023522 . 
90.85 51438826115 . - 30023529 . 
97.34 51438818018 . - 30023539 . 
103.83 51435956478 . - 30026876 . 
110.32 51434560951 . - 30028501 . 
116.81 51434837100 . - 30028179 . 
123.30 51435080714 . - 30027896 . 
129.79 51435311164 . - 30027628 . 
136.28 51435512393 . - 30027394 . 
142.76 51435699197 . - 30027178 . 
149.25 51435880508 . - 30026968 . 
155.74 51436047241 . - 30026776 . 
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25.00 


6654.639 


7.7294 


770.60 


26.00 


6655.596 


7.7283 


806.27 


27.00 


6655.751 


7 . 72 S 1 


841.63 


28.00 


& 655.108 


7.7289 


876.99 


29.00 


6653.711 


7.7306 


911.83 


30.00 


6651.615 


7.7330 


947.19 


31.00 


6648.943 


7.7362 


981.93 


32.00 


6645.784 


7.7399 


1016.57 


33.00 


6642.315 


7.7440 


1051.21 


34.00 


6638.654 


7.7483 


1085.74 


35.00 


6634.978 


7.7526 


1120.38 


36.00 


6631.441 


7.7568 


1154.81 


37.00 


6628.199 


7.7606 


1189.65 


38.00 


6625.379 


7.7639 


1223.36 


39.00 


6623.118 


7.7666 


1255.62 


40.00 


6621.515 


7.7684 


1286.84 


41.00 


6620.631 


7.7695 


1317.86 


42.00 


6620.505 


7.7696 


1348.26 


43.00 


6621.127 


7.7689 


1379.38 


44.00 


6622.477 


7.7673 


1410.30 


45.00 


6624.516 


7.7646 


1436.45 


46.00 


6627.376 


7.7609 


1462.20 


47.00 


6630.657 


7.7571 


1496.42 


48.00 


6634.186 


7.7530 


1529.72 


49.00 


6637.818 


7.7488 


1563.32 


50.00 


6641.411 


7.7446 


1597.65 


51.00 


6644.824 


7.7407 


1631.46 


52.00 


6647.925 


7.7371 


1664.65 


53.00 


6650.596 


7.7340 


1698.46 


54.00 


6652.744 


7.7315 


1732.27 


55.00 


6654.278 


7.7298 


1765.66 


56.00 


6655.151 


7.7288 


1799.68 


57.00 


6655.343 


7.7286 


1833.38 


58.00 


6654.841 


7.7292 


1867.40 


59.00 


6653.669 


7.7306 


1901.31 


60.00 


6651.866 


7.7327 


1934.71 


61.00 


6649.535 


7.7354 


1969.04 


62.00 


6646.736 


7.7387 


2003.05 


63.00 


6643.614 


7.7424 


2036.66 


64.00 


6640.256 


7.7464 


2070.47 


65.00 


6636.822 


7.7504 


2103.76 


66.00 


6633.438 


7.7544 


2136.74 


67.00 


6630.226 


7.7582 


2170.24 



162.23 51436208297 . - 30026590 . 
168.72 51436378759 . - 30026392 . 
175.21 51436530808 . - 30026217 . 
181.70 51436685225 . - 30026038 . 
188.19 51436850956 . - 30025846 . 
194.68 51437006143 . - 30025666 . 
201.17 51437196213 . - 30025446 . 
207.66 51437368920 . - 30025245 . 
214.15 51437587533 . - 30024991 . 
220.64 51437799569 . - 30024744 . 
227.13 51438037086 . - 30024467 . 
233.62 51438289123 . - 30024173 . 
240.10 51438574149 . - 30023840 . 
246.59 51438836529 . - 30023533 . 
253.08 51438840830 . - 30023530 . 
259.57 51438834154 . - 30023540 . 
266.06 51438836523 . - 30023538 . 
272.55 51438836153 . - 30023539 . 
279.04 51438840521 . - 30023533 . 
285.53 51438837255 . - 30023536 . 
292.02 51436551488 . - 30026202 . 
298.51 51434160812 . - 30028987 . 
305.00 51434443176 . - 30028656 . 
311.49 51434694318 . - 30028363 . 
317.98 51434926906 . - 30028091 . 
324.47 51435148627 . - 30027833 . 
330.95 51435348781 . - 30027600 . 
337.44 51435535867 . - 30027382 . 
343.93 51435712082 . - 30027177 . 
350.42 51435892956 . - 30026967 . 
356.91 51436042361 . - 30026794 . 
363.40 51436201010 . - 30026610 . 
369.89 51436371724 . - 30026411 . 
376.38 51436530981 . - 30026227 . 
382.87 51436693930 . - 30026037 . 
389.36 51436843665 . - 30025863 . 
395.85 51437033997 . - 30025642 . 
402.34 51437203180 . - 30025445 . 
408.83 51437417724 . - 30025195 . 
415.32 51437615550 . - 30024965 . 
421.80 51437843332 . - 30024699 . 
428.29 51438096137 . - 30024404 . 
434.78 51438361782 . - 30024093 . 
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68.00 


6627.315 


7.7616 


2203.95 


441.27 


51438659553 . 


- 30023745 . 


69.00 


6624.808 


“’.7646 


2236.21 


447.76 


51438839896 . 


- 30023534 . 


70.00 


6622 . S 32 


7.7669 


2267.43 


454.25 


51438841147 . 


- 30023534 . 


71.00 


6621.455 


7.7685 


2298.14 


460.74 


51438854059 . 


- 30023520 . 


72.00 


6620.730 


7.7694 


232 S .33 


467.23 


51438848302 . 


- 30023528 . 


73.00 


6620.676 


7.7694 


2358.32 


473.72 


51438838835 . 


- 30023539 . 


74.00 


6621.291 


7.7687 


2389.23 


480.21 


51438840540 . 


- 30023537 . 


75.00 


6622.551 


7.7672 


2419.84 


486.70 


51438841940 . 


- 30023534 . 


76.00 


6624.435 


7.7647 


2446.10 


493.19 


51436666305 . 


- 30026071 . 


77.00 


6627.090 


7.7612 


2471.84 


499.68 


51433850507 . 


- 30029352 . 


78.00 


6630.147 


7.7576 


2505.45 


506.17 


51434145633 . 


- 30029006 . 


79.00 


6633.448 


7.7538 


2538.84 


512.66 


51434396247 . 


- 30028713 . 


80.00 


6636.870 


7.7499 


2572.24 


519.14 


51434642574 . 


- 30028425 . 


81.00 


6640.288 


7.7459 


2605.53 


525.63 


51434851838 . 


- 30028181 . 


82.00 


6643.582 


7.7421 


2639.14 


532.12 


51435066093 . 


- 30027931 . 


83.00 


6646.639 


7.7386 


2672.94 


538.61 


51435260419 . 


- 30027705 . 


84.00 


6649.350 


7.7354 


2706.03 


545.10 


51435436862 . 


- 30027500 . 


85.00 


6651.626 


7.7328 


2739.84 


551.59 


51435610960 . 


- 30027297 . 


86.00 


6653.393 


7.7308 


2773.55 


558.08 


51435778411 . 


- 30027103 . 


87.00 


6654.598 


7.7294 


2806.74 


564.57 


51435948364 . 


- 30026905 . 


88.00 


6655.195 


7.7287 


2840.13 


571.06 


51436102292 . 


- 30026727 . 


89.00 


6655.170 


7.7288 


2874.05 


577.55 


51436255865 . 


- 30026549 . 


90.00 


6654.533 


7.7295 


2907.24 


584.04 


51436418385 . 


- 30026360 . 


91.00 


6653.305 


7.7310 


2940.53 


590.53 


51436578220 . 


- 30026174 . 


92.00 


6651.546 


7.7331 


2974.03 


597.02 


51436757933 . 


- 30025965 . 


93.00 


6649.293 


7.7357 


3007.74 


603.51 


51436918736 . 


- 30025778 . 


94.00 


6646.656 


7.7388 


3041.55 


609.99 


51437110521 . 


- 30025555 . 


95.00 


6643.718 


7.7423 


3075.35 


616.48 


51437314954 . 


- 30025317 . 


96.00 


6640.575 


7.7460 


3109.16 


622.97 


51437527062 . 


- 30025070 . 


97.00 


6637.340 


7.7498 


3142.25 


629.46 


51437748538 . 


- 30024811 . 


98.00 


6634.132 


7.7536 


3175.44 


635.95 


51437987351 . 


- 30024532 . 


99.00 


6631.065 


7.7572 


3208.84 


642.44 


51438257984 . 


- 30024216 . 


100.00 


6628.231 


7.7605 


3242.54 


648.93 


51438533304 . 


- 30023894 . 
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